#include <math.h>
#include <R.h>
#include "geom3.h"
#include "functable.h"

/*
	$Revision: 1.1 $	$Date: 2009/11/04 23:54:15 $

	K function of 3D point pattern


	k3trans	  	translation correction

	k3isot		isotropic correction


# /////////////////////////////////////////////
# AUTHOR: Adrian Baddeley, CWI, Amsterdam, 1991.
#         
# MODIFIED BY: Adrian Baddeley, Perth 2009.
#
# This software is distributed free
# under the conditions that
# 	(1) it shall not be incorporated
# 	in software that is subsequently sold
# 	(2) the authorship of the software shall
# 	be acknowledged in any publication that 
# 	uses results generated by the software
# 	(3) this notice shall remain in place
# 	in each file.
# //////////////////////////////////////////////


*/


void
k3trans(p, n, b, k)
     Point *p;
     int n;
     Box *b;
     Ftable *k;
{
  register int i, j, l, lmin;
  register double dx, dy, dz, dist;
  register double  vx, vy, vz;
  Point *ip, *jp;
  double dt, vol, lambda, denom, term;
  double sphesfrac(), sphevol();

  /* compute denominator & initialise numerator*/
  vol = (b->x1 - b->x0) * (b->y1 - b->y0) * (b->z1 - b->z0);
  lambda = ((double) n )/ vol;
  denom = lambda * lambda;

  for(l = 0; l < k->n; l++) {
    (k->denom)[l] = denom;
    (k->num)[l]   = 0.0;
  }

  /* spacing of argument in result vector k */
  dt = (k->t1 - k->t0)/(k->n - 1);

  /* compute numerator */
  for( i = 0; i < n; i++) {
    ip = p + i;
    for(j = i + 1; j < n; j++) {
      jp = p + j;
      dx = jp->x - ip->x;
      dy = jp->y - ip->y;
      dz = jp->z - ip->z;
      dist = sqrt(dx * dx + dy * dy + dz * dz);
      lmin = ceil( (dist - k->t0) / dt );
      if(lmin < 0)
	lmin = 0;
      
      vx = b->x1 - b->x0 - (dx > 0 ? dx : -dx);
      vy = b->y1 - b->y0 - (dy > 0 ? dy : -dy);
      vz = b->z1 - b->z0 - (dz > 0 ? dz : -dz);
      if(vx >= 0.0 && vy >= 0.0 && vz >= 0.0) {
	  term = 2.0 /(vx * vy * vz);
	  /* 2 because they're ordered pairs */
	  for(l = lmin; l < k->n; l++)
	    (k->num)[l] += term;
      }
    }
  }

  /* compute ratio */
  for(l = 0; l < k->n; l++)
    (k->f)[l] = ((k->denom)[l] > 0.0)?
      (k->num)[l] / (k->denom)[l]
      : 0.0;
}

void
k3isot(p, n, b, k)
     Point *p;
     int n;
     Box *b;
     Ftable *k;
{
  register int i, j, l, lmin;
  register double dx, dy, dz, dist;
  Point *ip, *jp;
  double dt, vol, denom, term;
  double sphesfrac(), sphevol();
  Point vertex;
  Box   half;

  /* compute denominator & initialise numerator*/
  vol = (b->x1 - b->x0) * (b->y1 - b->y0) * (b->z1 - b->z0);
  denom = ((double) (n * n))/vol;

  for(l = 0; l < k->n; l++) {
    (k->denom)[l] = denom;
    (k->num)[l]   = 0.0;
  }

  /* spacing of argument in result vector k */
  dt = (k->t1 - k->t0)/(k->n - 1);

  /* set up for volume correction */

  vertex.x = b->x0;
  vertex.y = b->y0;
  vertex.z = b->z0;
  half.x1  = b->x1;
  half.y1  = b->y1;
  half.z1  = b->z1;
  half.x0  = (b->x0 + b->x1)/2.0;
  half.y0  = (b->y0 + b->y1)/2.0;
  half.z0  = (b->z0 + b->z1)/2.0;

	/* compute numerator */
  for( i = 0; i < n; i++) {
    ip = p + i;
    for(j = i + 1; j < n; j++) {
      jp = p + j;
      dx = jp->x - ip->x;
      dy = jp->y - ip->y;
      dz = jp->z - ip->z;
      dist = sqrt(dx * dx + dy * dy + dz * dz);
      lmin = ceil( (dist - k->t0) / dt );
      if(lmin < 0)
	lmin = 0;
      
      term = (1.0 / sphesfrac(ip, b, dist)) 
	+ (1.0 / sphesfrac(jp, b, dist)); 
      term *= 
	1.0 - 8.0 * sphevol(&vertex, &half, dist) / vol;
      for(l = lmin; l < k->n; l++)
	(k->num)[l] += term;
    }
  }
  /* compute ratio */
  for(l = 0; l < k->n; l++)
    (k->f)[l] = ((k->denom)[l] > 0.0)?
      (k->num)[l] / (k->denom)[l]
      : 0.0;
}
